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We study the initial instability of flat sand surface and further nonlinear dynamics of wind 
ripples. The proposed continuous model of ripple formation allowed us to simulate the development 
of a typical asymmetric ripple shape and the evolution of sand ripple pattern. We suggest that this 
evolution occurs via ripple merger preceded by several soliton-like interaction of ripples. 



I. INTRODUCTION 

Aeolian (i.e. wind driven) sand ripples form nice regular patterns on coastal beaches and desert floors and indicate 
an instability of flat sand surface under the wind-induced transport and rearrangement of loosely packed sand grains. 
Following fundamental work by Bagnold § more than half a century ago, formation of sand ripples has been studied 
by many researchers (see review Q and the references therein). Significant progress in understanding the nature of 
t~* this phenomenon has been achieved. Nonetheless, major questions remain open; these involve the most interesting 
part of ripple formation, the nonlinear interactions that follow the initial instability. Previous research on ripples has 
generally relied on highly simplified continuum models or on stochastic or molecular dynamics simulation. By means 
of a new deterministic continuous model that seems to better describe the essential physics, we here investigate salient 
nonlinear properties of ripple formation. 

As is well known, Aeolian ripples are oriented perpendicularly to the wind direction. Mature ripples are asymmet- 
rical in cross section: their stoss (upwind) slopes are typically much less steep than the shorter lee (downwind) slopes 
The steepness of the lee slopes cannot exceed and usually does not reach the sand angle of repose. The ripples 
have convex stoss slopes, concave lee slopes, and flattened crests which usually end with a brink. Since smaller ripples 
of the same shape have smaller volume to surface ratios, they are translated faster by the wind and can overtake the 
larger ripples. A possible merger results in gradual elimination of small ripples and in growth of ripple wavelength. 

To analyze the mechanics of sand transport, which occurs whenever the wind is sufficiently strong, it is convenient 
to distinguish two types of sand grain movement: saltation and reptation (or creep) Jl|,||. Saltating grains move 
by long trajectories that end in high-energy impacts with the surface. These impacts take place at almost uniform 
• shallow angles of descent varying from 9 to f5° to the horizontal depending on the wind strength and grain size [|J. 
After an impact, a saltating grain usually rebounds sufficiently high to be accelerated by the wind again and continues 
its saltation. These grains gain energy from the wind and transfer part of it to the sand bed on impact. Each impact 
may cause ejection of one or several low energy (reptating) particles from the bed surface. Reptating particles make 
a short hop, usually jumping or rolling for several millimeters or less, and stop. The exchange flux between saltating 
and reptating grain populations is supposed to be small HQ . 

According to the hypothesis put forward by Bagnold [Q , the ripple wavelength is equal to the mean length of saltation 
jump. However, this claim has been challenged by numerous researchers (see, e.g., [|]](|) and it is now commonly 
accepted that the essential physics lies in the variation of reptation flux. The role of saltation, whose trajectories 
•i-H . are many times longer than the ripple wavelength, is indirect. The oblique, almost unidirectional bombardment by 
' saltating particles supplies the energy necessary for reptation. Wind strength determines the intensity of saltation; 
the probability of direct entrainment of particles into reptation by wind is small Q| . 

A model, based on these views, was proposed by Anderson ||. Linear stability analysis of this model showed 
that the initial ripple wavelength is determined by, and is several times larger than, the mean length of reptation. 
Unfortunately, the model yields unrealistic results at the nonlinear stage of ripple growth, which begins very early. 
It has been suggested Q that the model can be improved by allowing the reptating grains to continue rolling upon 
the bed surface after landing, and not to stop immediately as was assumed originally. Although no such continuous 
model has been developed, employing a similar approach in molecular dynamics computer simulations of sand ripples 
JjJ was quite successful. 

Interesting results on modeling different aspects of sand ripples dynamics have been obtained by means of cellular 
automata models or molecular dynamics simulations |7-10|. It was even claimed ]Io| that no continuum mechanics 
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or deterministic model can capture the main feature of ripple self-organization: the increase of scale in time due 
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to merging of ripples. Below, we show that this general conclusion, based on a schematic discrete model of ripple 
merging, is wrong. The continuous deterministic model proposed herein provides a better description of the physical 
process than the simplified stochastic model |Q . Our model reproduces the asymmetrical ripple shape and is able to 
simulate not only merging of ripples but also a more complicated, soliton-like mode of ripple interaction which can be 
equally significant for ripple self-organization. A somewhat similar behavior has been observed in molecular dynamics 
simulation of sand ripples fjj . 

Let us mention here also the analytical model of sand ripples . In this model, the upwind slope of a ripple is 
given by a smooth solution of the diffusion equation with a negative diffusion coefficient. Such a solution is, however, 
unstable. Furthermore, in the absence of wind the ripple shape in jnj is described by a well-posed diffusion equation. 
This is also unphysical: according to this model ripples diffuse and disappear, since surface particles continue to roll 
down the slopes. Real sand ripples are, of course, metastable. 

To derive a model allowing for metastability, it is necessary to take into account that the surface flux of granular 
material is not determined solely by the local surface slope. Models involving an additional variable, surface flux 
or the density of rolling particles, have recently been derived to simulate quasistationary evolution of a growing pile 
shap e jl2|l or to model the dynamics of pile surface in more detail and on a shorter spatio-temporal scale ( see 
also ]14[). In our model of Aeolian ripples we use a similar approach to describe the flow of rolling particles. 



II. MATHEMATICAL MODEL 



In accordance with the physical picture of sand transport described above, we assume that there exists a uniform 
flux of saltating grains which reach the bed surface y = /i(x, t) at a low angle 7 to horizontal. The impacts cause 
erosion of this surface. The erosion rate / is proportional to the impact intensity, which depends on the surface 
orientation with respect to the direction of saltation. Let the saltating particles strike an inclined surface at an angle 
0. Then / is proportional to smO and we can write 

sin# 
J — jo — 
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where /o is the rate of erosion of a horizontal surface (determined by the intensity of saltation). However, if the 
surface has sufficiently steep slopes, some parts of it may be in shadow and unreachable by saltating grains. In this 
case we set / = 0; shadowing introduces non-locality into this problem. 

A reptating particle, ejected by an impact at a point x, makes a jump and lands on the bed surface at a point y 
with probability density p given by the " splash function" , p = p a (x, y) , first introduced in [jl5| . This function, which 
will be specified below, depends, in particular, on the surface slope at the point of impact. In our model, the splash 
function also accounts for all the anisotropy induced by a chosen wind (saltation) direction. 

Upon landing, reptating particles do not stop immediately but may roll away, although usually not far from the 
landing point. Let i?(x,i) be the effective surface density of rolling particles (Rdx is the volume which particles, 
presently rolling over the part of the free surface above the area <ix, would occupy in the sand bed). When they 
stop, the rolling particles are incorporated into the motionless bed. Following p3| , we denote by T[h, R] the rate of 
rolling-to-steady state transition and write the mass conservation equations for the sand bed and for the population 
of rolling particles: 

d t h = T[h,R]-f, (1) 
d t R + V • J = Q - T[h, R]. (2) 

Here J is the horizontal projection of the flux of rolling particles, and the source term 

Q(x,i) = J f(y,t) Pa (y,x)dy (3) 

gives the intensity of "rain" of falling reptating particles. 

We assume that reptating particles lose most of their momentum in collision with the rough bed surface. Neglecting 
inertia, we postulate that upon landing the particles roll in direction of the steepest descent and that the steeper the 
slope the faster they roll. The simplest form [jl6| of the flux J is, therefore, 

J = —fioRVh, 

where fio is a constant "mobility" of particles. 
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The rate of rolling-to-steady state transition, T, depends on stability of a particle on inclined sand bed surface and 
on the amount of rolling particles. Rolling particles never form a thick layer on the surface during the ripple growth: 
there is only a small amount of them at each time moment. It can be assumed that, for a fixed free surface incline, 
the rate of rolling-to-steady state transition is proportional to the amount of rolling particles on the surface, R. Since 
the exchange rate cannot depend on the free surface slope orientation, we further assume T is a (smooth) function of 
|V h\ 2 . The steeper the free surface is, the easier do particles roll down and, correspondingly, the lower is the rate of 
rolling-to-steady state transition. For slopes steeper than the sand angle of repose, a r ~ 30°, rolling sand grains do 
not stop at all. Taking these arguments into account we assume, as a simple but physically reasonable approximation, 
that r is also proportional to tan 2 a r — \ V h\ 2 and obtain 



T = k R 1 



|Vfr| 2 
tan 2 a r 



where Ko characterizes particle stability upon a horizontal surface. 
We now rescale the variables, 



t' = f t, R' = ^R, /' = !/, J' = |j, r' = |r, 

Jo Jo Jo Jo 



and obtain 



J = -vRVh, r = R I 1 - -M- . . 

tan a r J 



(4) 



and / = sin 8/ sm^, or / = in a shadow. Here 9 — 9(x,t) is the angle at which the saltating particles strike the 
bed, and v = (j-o/kq is a constitutive dimensionless parameter characterizing the competition between mobility and 
stability of dislodged grains on the sand surface. The rescaling leaves the equation (|l) invariant while the equation 
(0) takes the form 



^-dtR + V -J = Q-T[h,R\. 

K 



(5) 



We suppose that the rate of sand surface erosion caused by saltation, /o, is usually much smaller than kq, the 
coefficient determining the rate at which rolling particles come to a stop on the horizontal surface and are absorbed 
by the motionless bulk (that is why the layer of rolling particles is so thin) . Neglecting the small term in equation (|^) 
we arrive at a quasistationary mass balance equation for rolling particles, 



V ■ J = Q-T[h,R]. 



(G) 



III. SPLASH FUNCTION 

To complete the model (Q), (||), ([I]), and (|^) we need to specify the splash function. Although not much is known 
about this function, previous studies (see, e.g., ||-[8]]) suggest that the system is not very sensitive to the details of 
splash function behavior and that an approximation, sufficient at least for qualitative simulation, may be obtained 
by combining the existing experimental data and simple physical arguments. We limit our consideration to the 
one-dimensional (ID) case. 

Collision of quartz grains with a sand bed has been studied experimentally by Willetts and Rice ||. It was found 
that ejection of reptating grains from the bed depends only slightly on the incident angle of attack, which varied in 
these experiments in accordance with the systematic changes of the bed inclination. (The angle of descent of saltating 
particles to the horizontal was constant). For various incident angles, ejection occurred at approximately the same 
mean angle to the bed surface, w 50°, with standard deviation w 40°, and the same mean velocity of ejecta. 

We use these results to crudely reconstruct the dependence of the splash function on the bed surface inclination 
in ID case. First, we define (somewhat arbitrarily, see fl%) the density of ejection angle distribution, p = p((j>), 
providing for the mean and standard deviation values as found in |^| (see Fig. la). For simplicity, we further assume 
that particles are ejected from the bed with the same initial velocity vq at different angles 4>, that they then follow 
simple ballistic trajectories x = xq + vq cos(a + (j))t, y — h(xo) + VQ sin(a + 0)£— ^gt 2 j an d hit the surface at a horizontal 
distance s from the ejection point. Here a = t&n~ 1 (d x h) is the surface angle at the ejection point xq, g - acceleration 
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of gravity. Assuming the surface curvature is small, |<3| x /i| <C q/vq, we approximate h(x$ + s) by h(xo) + sd x h(xo) 
and find the jump length 

v 2 

s = —(sin 2{a + (f\ — 2 cos 2 {a + </>} tana). 

Making use of the probability density p{4>) , we can now calculate numerically the mean and the standard deviation 
of reptation jump length, m r (a) and a r (a), for any bed surface inclination a, up to the value of a proportionality 
coefficient v^/g. This factor is eliminated from the final dimensionless formulation of the model by choosing the unit 
of length equal to the mean reptation length at the horizontal bed surface. Thus m r (0) = 1 by definition. The 
functions m r (a) and ay (a) are shown in Fig. fb. Finally, we approximate the splash function p a by the density of a 
corresponding normal distribution, 
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where a = a r (a), m = m r (a), and a — tan 1 d x h(xo,t). 



IV. LINEAR STABILITY ANALYSIS 



To analyze the initial instability of a horizontal sand bed, we assume that h and its derivatives are small and 
linearize the model. Up to the second order terms, T[h, R] — R and / = f + kjd x h, where fc 7 = cot 7, so Eq. (|l|) 
yields R = 1 + d t h + k 7 d x h. 

For small surface slopes, the standard deviation of reptation length does not change much (see Fig. lb) and we 
set a r = oy(0) = 1.25. Knowing the dependence m r (a) (Fig. lb), it is easy to find numerically that for small slopes 
m r 1 — k m d x h, where k m = 2.01. 

Let po(x) be the density of the normal distribution iV(l,a 2 (0)). It is not difficult to show that Pa{x,y) = pa(y — 
x) + k m p' (y — x)d x h{x) plus the higher order terms (""' means derivative). The linearized Eq. (|3|) takes the form 
Q = 1 + kypo * d x h + k m p' * d x h, where "*" is the operator of convolution. 

Substituting the linear approximations for T, R, and Q into Eq. (^) we obtain, up to the second order terms, 

d t h = vd 2 xx h + fc 7 (p * d x h - d x h) + k m p' * d x h 

We can now apply the Fourier transform and substitute h = e xt+lux to find the dispersion relation 

A(o>) = —vu 2 + kjiw(po — 1) — k m ui po, 

where po — exp(— [wa r (0)] 2 /2 — iw) is the Fourier transform of po. Note that with v = k m = one gets the dispersion 
relation for Anderson's model fl. The initial ripple wavelength can be calculated as Iq = 2it/ujq, where ujq is the wave 
number at which the expression 

Re\{ijj) = -vuj 2 + w(fe 7 sinw - k m uj cos w )e~ (t ^ (0))2/2 

attains a positive maximum. This maximum exists and the flat surface is unstable if 

kj > k m + v. (7) 

To explain this result we note that ripples grow because of the geometrical effect of greater impact and ejection flux 
on upwind-oriented slopes than on downwind-oriented slopes Q . This non-uniformity increases if the saltation angle 
of attack becomes smaller (fc 7 increases). On the other hand, the greater v and k m are, the more significant are, 
respectively, smoothing effects due to rolling of dislodged particles down the surface slopes and scattering of ejected 
particle trajectories by an uneven sand surface. Condition (^) indicates when the instability prevails. 

As it was mentioned above, the saltating particles usually hit the sand surface at almost uniform angle to the 
horizontal, varying from 9 to 15°. In the examples below we use the mean reported value of this angle, 7 = 12°. 
Using the stability condition (0), it is easy to calculate that the flat surface is unstable if v < 2.69. 
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V. NONLINEAR DYNAMICS 



To study the ripple evolution further we solved the nonlinear system (Q), |^), (Q), (^) numerically, assuming periodic 
boundary conditions and using an implicit finite-difference approximation. The initial evolution of a slightly disturbed 
flat surface obeys the linear theory: after a short initial stage the fastest growing mode, having the wavelength 
predicted by the linear stability analysis, dominates. 

The growing ripples remain almost symmetric until, shortly before the appearance of first shadow regions, the 
downwind slopes become steeper. This asymmetry develops quickly as the ripples continue to grow. Further evolution 
is accompanied by the increase of the ripple wavelength (see Fig. 2). It can be seen that the downwind translation 
of ripples gradually slows down as their size increases. Although our model is much simplified in many respects, the 
calculated mature ripple shape is similar to that of the real sand ripples for v = 2 wc obtained convex stoss 
slopes inclined at about 11 — 13°, flattened crests, and slightly concave lee slopes with the mean maximal inclination 
26 — 27°; the ripple index (length to height ratio) was about 14. For v = 1 the results are qualitatively similar, 
although the mean maximal inclinations are 15° and 28-29° for the stoss and lee slopes, respectively, and the ripple 
index is smaller. 

The most interesting part of ripple dynamics is the mechanism of ripple merging and self-organization. The simu- 
lations show that simple merger takes place only if the overtaking ripple is much smaller than the bigger " overtaken" 
ripple, which moves more slowly. Otherwise another, more complicated scenario is usually realized. As a smaller 
ripple reaches the larger, the trough between them becomes shallow and a "two-headed" long ripple appears. Then 
the "downwind head", which originally formed the larger ripple crest, starts to move forward as a separate small 
ripple and runs ahead. Two new ripples emerge from this recombination. The ripple that is left behind is larger than 
the larger of the two ripples before merger. The ripple that runs away is smaller than the smaller one before this event 
(Fig. 3). Sometimes the trough between the two merging ripples disappears before the generated long ripple becomes 
unstable. However, soon there appears a new trough near the end of this ripple crest and a small running away 
ripple develops. Since it is smaller, this ripple proceeds even faster and soon meets another large ripple on its way. 
Complete merger is now more probable, since the overtaking ripple became smaller. However, another recombination 
may yet occur before the material redistribution between ripples is completed. Obviously, such a mechanism of ripple 
interaction produces a regular ripple array structure more efficiently than a simple merger of ripples. 

In our opinion, this scheme of ripple interaction gives also a likely explanation to the appearance of small secondary 
ripples in wind tunnel experiment |0|. Indeed, the appearance of such ripples due to the backward eddy flow behind 
a ripple, as is suggested in jlfifl , seems hardly possible. Sand ripples are so shallow that most probably there is no 
backward flows in their shadows. Even if small backward eddies existed, they could never cause saltation of sand 
grains against the main wind direction, and thus could not produce any sand ripples. As follows from our simulation, 
no backward flows are necessary: small ripples appear as a means of redistribution of mass during the ripple array 
reconstruction leading to the wavelength growth. 
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FIG. 1. a - density of ejection angle distribution, p(<f>); b - the mean (solid line) and standard deviation (dashed line) of 
reptation jump length as functions of surface inclination. 

FIG. 2. Flat surface instability and formation of sand ripples. Wind direction is from left to right; the unit of length is the 
mean length of the reptation jump; v — 2. 

FIG. 3. A typical ripple interaction; see the region bounded by the thin line. To show the details the ripples are stretched 
in the vertical direction. 
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